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Two aspects of TDDFT, the linear response approach and the adiabatic local density approx- 
imation, are examined from the perspective of lattice models. To this end, we review the DFT 
formulations on the lattice and give a concise presentation of the time-dependent Kadanoff-Baym 
equations, used to asses the limitations of the adiabatic approximation in TDDFT. We present 
results for the density response function of the 3D homogeneous Hubbard model, and point out a 
drawback of the linear response scheme based on the linearized Sham-Schliiter equation. We then 
suggest a prescription on how to amend it. Finally, we analyze the time evolution of the density in a 
small cubic cluster, and compare exact, adiabatic-TDDFT and Kadanoff-Baym-Equations densities. 
Our results show that non-perturbative (in the interaction) adiabatic potentials can perform quite 
well for slow perturbations but that, for faster external fields, memory effects, as already present in 
simple many-body approximations, are clearly required. 
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I. SCOPE OF THIS WORK 

A detailed knowledge and manipulation of non- 
equilibrium phenomena is one of the great challenges of 
current research in condensed matter - both from a fun- 
damental point of view and from the perspective of pos- 
sible applications. Theoretical research can contribute 
significantly to this endeavor, by answering a number of 
important questions. This is what motivates the strong 
effort in developing theoretical methods to describe sys- 
tems out of equilibrium. 

In this work we focus on one of these methods, namely 
time-dependent density- functional theory (TDDFT) [T]. 
TDDFT provides the extension to the time-dependent 
case of static density-functional theory (DFT) [21 13] . Its 
conceptual foundations were laid in the mid-eighties with 
the Runge-Gross theorem [3]. Since then, alternative 
proofs of such theorem have been presented [5] U (see 
also ^), together with variants for ensembles (8j, multi- 
component systems [51, open systems [TUl, superconduc- 
tivity [TTl [T^l , or for the case of time-dependent current- 
density-functional theory (TDCDFT) [13-17 ■ Further- 
more, connections with other non-equilibrium methods 
have also been examined to clarify general conceptual 
issues in TDDFT ^\T^. 

In TDDFT, the key variable is the one-particle den- 
sity n, and a central ingredient is the time-dependent 
exchange-correlation (XC) potential Vxc^ which incorpo- 
rates all the complexities of the time-dependent many- 
body dynamics. As a result of this contracted descrip- 
tion, and since time enters explicitly into the formulation, 
Vxc depends in a highly non-trivial way on the entire his- 
tory of the density n (memory effects). TDDFT is, in 
principle, an exact scheme. In practice, a major disad- 
vantage in concrete applications is the lack of an accurate 
description of dynamical inter-particle correlations. 

An easy but inadequate way to proceed is to make 
use of the so-called adiabatic local density approximation 



( ALDA) f5D] , where the XC potential at every particular 
time only depends on the local density. This amounts 
to neglecting non-locality in space and memory effects in 
Vxc- In general, in order to improve the ALDA, it be- 
comes necessary to consider the issue of ultranonlocality 
[?T1E51 . This means that the introduction of non-locality 
in time requires the inclusion of strong non-local effects 
in space. 

Some simplifications can be made in the case of weak 
perturbing fields, where linear response arguments apply 
[21] . However, also in this case, the inclusion of non-local 
effects is a far-from-trivial task, as further discussed in 
the rest of the paper. Being a direct gateway to the 
study of important properties such as, e.g., the optical 
response of materials, linear response within TDDFT has 
received a great deal of attention |25]. The same applies 
to the problem of improving on the ALDA away from the 
linear regime: there has been a considerable theoretical 
effort in developing reliable XC potentials for far-from- 
equilibrium dynamics [26H37j . where memory effects are 
important. 

In this work we analyze these two issues from a quite 
specific perspective, i.e. by considering TDDFT applied 
to lattice model systems. Lattice models have a long his- 
tory in condensed matter physics, due to their conceptual 
simplicity and to the fact that many of such systems of- 
ten admit an exact analytic solution. It thus does not 
come as a surprise that lattice models have also been in- 
vestigated in the context of DFT and TDDFT in order 
to provide theoretical insight into fundamental aspects of 
these frameworks. 

The plan of our paper is as follows. We start in Sect. |ll] 
with a presentation of the linear response formalism, fol- 
lowed by a review of (TD)DFT for lattice systems in Sect. 



Ill This, especially for the time-dependent case 



is a rather new topic, and we will also survey the recent 
literature in the field. Then, in Sect. |IV| we proceed 
to a short summary of another method for treating non- 
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equilibrium problems - the solving of the KadanofF-Baym 
equations (KBE) |40l [41] . Here approximations for elec- 
tron correlations are constructed from standard meth- 
ods of many-body perturbation theory (MBPT). Results 
from this approach will be used as benchmark in our 
analysis of the ALDA for lattice models, which is done 
in Sect. |Vj In this section, we will also illustrate, using 
small model linear chains, some of the general aspects 
of linear response theory within TDDFT, as discussed in 
Sect, [n] Finally, our conclusions are provided in Sect. 

ED 



II. LINEAR RESPONSE 

First-order perturbation theory and, what amounts to 
something very similar in spirit, the theory of linear re- 
sponse has been around since the advent of quantum me- 
chanics. Unfortunately, first-order perturbation theory is 
rarely sufficient and more and more complicated higher 
order corrections must be included. In infinite systems 
summations to infinite order of carefully selected correc- 
tions must be carried out. Perturbation expansions do, 
in principle, not converge. They are, at best, asymptotic 
in character and physical intuition reigns - even in finite 
systems as has been realized in recent years [42]. Yet, 
perturbation theory often provides valuable insight into 
the qualitative deviations between the behavior of sim- 
ple models and the realistic systems they are intended to 
describe. 

Another very important use of perturbation theory is 
the description of weak experimental probes applied to 
systems for the purpose of investigating their properties 
without seriously affecting those very properties. A typi- 
cal example is the forced interaction between matter and 
weak electromagnetic radiation. If the intensity of the 
radiation is chosen weak enough, first-order perturbation 
theory suffices to map out the experimental results which 
are then described in terms of so called linear response 
functions like, e.g.^ the density-density or the current- 
current response functions. 

In this section we will be mainly concerned with the 
density-density response function which describes the 
first-order response of the system to a perturbing local 
potential. The excitation of particle-hole pairs in atoms, 
molecules, and solids as well as the absorption of light by 
such systems are well within the realm of this theory. 

The traditional way of constructing the density re- 
sponse function of an electronic system is to include the 
long-range Coulomb interaction through infinite-order 
perturbation theory. A very successful way of accounting 
for the important interaction between the excited elec- 
tron and the associated hole it leaves behind is the Bethe- 
Salpeter method [HJSnj. Here, the vertex corrections are 
obtained by solving a two-particle Schrodinger-like equa- 
tion for the particle and the hole interacting through a 
statically screened Coulomb interaction. This method 
has been very successful in describing the strong excitonic 



effects in, e.g., rare-gas systems [H]. Unfortunately, the 
method is computationaly very demanding and the treat- 
ment of systems with low symmetry like surfaces or large 
molecules are often beyond reach. 

In the last decade, time-dependent density-functional 
theory |3J H?) has emerged as a tool for calculating the 
density responses of realistic systems. The great compu- 
tational advantage of TDDFT is that the vertex correc- 
tions are much simpler two-point correlation functions. 
The drawback, however, is the lack of a systematic ap- 
proach to finding successively better approximations to 
the exchange-correlation kernel which contains all the 
important particle-hole interactions beyond the random 
phase approximation (RPA). Within TDDFT the full 
density response function x is given by |24j 

X^Xo + Xo{v + fxc)x (1) 

where xo is the non-interacting density response function 
expressed in Kohn-Sham one-electron orbitals, v is the 
bare Coulomb interaction (1/r), and fxc is the famous 
exchange-correlation kernel into which all the beyond- 
RPA many-body effects have been deferred. The ker- 
nel fxc is formally defined to be the second functional 
derivative of the XC part of the action functional of 
TDDFT. A common formal starting point for obtain- 
ing approximations to this kernel is the so called Sham- 
Schliiter (SS) [15] equation expressing the fact that the 
densities of the fully interacting system and the Kohn- 
Sham system must be the same. If G is the exact one- 
electron Green's function of the interacting system and 
Go is the non-interacting Green's function of the equiva- 
lent Kohn-Sham system these Green's functions are con- 
nected through a Dyson-like equation 

G = Go + Go(E[G]-^;,,)G (2) 

where S is the self-energy of the interacting system and 
Vxc is the XC part of the local one-electron potential 
which generates Gq. Since the particle density of the 
system is given by the trace of the Green's function, tak- 
ing the trace of both sides of the Dyson's equation gives 
the SS: 

Tr[GGo«,c] =Tr[GI][G]Go] (3) 

The SS is thus an equation for the determination of the 
exact Vxc of TDDFT. A further variation with respect 
to, e.g., the external potential provides an expression for 
the XC kernel fxc = 5vxc/^n although it is rather com- 
plicated and has so far, to our knowledge, not been used 
by anyone for practical calculations. 

By doing ordinary diagramatic perturbation expan- 
sions for the Green's function G and the corresponding 
self-energy E in the SS, this equation provides a natural 
connection between TDDFT and many-body perturba- 
tion theory. Along these lines, a diagrammatic procedure 
for finding better kernels fxc'-^ was developed in Ref. [43j . 
but it still remains untested. 
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What has been used in several applications is a simpli- 
fied version of the SS called the linearized Sham-Schliiter 
equation (LSS). This is obtained by replacing every inter- 
acting Green's function G by the non-interacting Kohn- 
Sham Green's function Go in the SS. The resulting equa- 
tion 



Tr[GoGot;.c] = Tr[GoS[Go]Go] 



(4) 



can also be derived from the variational many-body ap- 
proach |35j. One then starts from the so called Klein 
[50j functional for the total action. Subsequently, one 
chooses to restrict the normally free variations of the 
one-electron Green's function G to those non-interacting 
ones that can be generated by multiplicative, local one- 
body potentials. This demonstrates that the LSS is a 
more accurate expression for generating XC potentials of 
TDDFT than is suggested by its original derivation where 
it was simply the lowest order result in a perturbation 
expansion. We notice that, within the LSS approach, 
very high-order correlation effects can be accounted for 
by choosing a very sophisticated self-energy. We should, 
however, also point out that there is an inherent lack of 
self-consistency in the LSS approach. The Kohn-Sham 
Green's function Go is certainly not generated by that 
Dyson equation which has S[Go] as a self-energy. This 
lack of self-consistency could easily lead to violations of 
some of the conservation laws and consistency require- 
ments which one usually takes for granted within varia- 
tional formulations. 

Unfortunately, there is a much more severe problem 
associated with use of the LSS approach as we will now 
discuss. One can easily perform a variation with respect 
to the externally applied potential in the LSS. Contain- 
ing, however, only the non-interacting Green's function 
Go this is equivalent to performing a variation 5V in the 
full effective Kohn-Sham potential V. Using the obvious 
fact that 



(5Go = GqSVGq 



(5) 



such a variation leads to the following equation for the 
XC kernel f^c 



d2d2' xo(l,2)/.c(2,2')xo(2',l') = 
i I d2dMAdb Go(l,2)Go(3,l)|^7^Go(l',5)Go(4,l') 



(6) 



where the quantity A is the difference between the chosen 
self-energy E and the XC potential v^c 



\5G(4,5) 

I I d2d2' Go(l,l')Go(l',2)A(2,2')Go(2',l) 



I / d2d2' Go(l,2)A(2,2')Go(2',l')Go(l',l) 



A(l,l') = S(l,l')-z;.c(l)'5(l,l') 



(7) 



If the four-point vertex function ST,/6G is taken from 
Hartree-Fock theory, which means that it becomes the 



bare Coulomb interaction, the resulting expression for 
the XC kernel fxc is called the exact-exhange approx- 
imation (EXX) within TDDFT. In the course of time, 
this approximation has been derived by many researchers 
and it has been given many different names like the opti- 
mized potential method (0PM), the optimized effective 
potential (OEP) approach, or the exchange-only approx- 
imation (EOA). It has been shown to have many nice 
properties in extended systems and it gives very accurate 
total energies in both finite and infinite systems ranging 
from atoms and molecules to the electron gas. The spec- 
tral properties of the resulting density response function 
was rather recently investigated at length in series of pa- 
pers [51H53]. 

In these works, it was discovered that the resulting re- 
sponse function of finite systems actually has poles in the 
upper half plane thus precluding a resonable description 
of the optical response of the system at higher energies. 
This problem is associated with rather unexpected ze- 
roes of the non-interacting response function xo st cer- 
tain frequencies. The presence of such zeroes was, how- 
ever, established a long time ago j53]. From the defining 
equation above for the kernel fxc we see that one has 
to invert the response function xo twice in order to ob- 
tain fxc- Thus, a zero in xo produces a double pole in 
fxc which, in turn, causes the full response function x 
to have an unphysical pole in the upper half plane. It 
is important to notice that this problem has nothing to 
do with the degree of sophistication by which one tries 
to incorporate the correlation effects. Any choice of self- 
energy will produce the same unphysical result. 



III. DFT AND TDDFT FOR LATTICE MODELS: 
THE CASE OF THE HUBBARD MODEL 

The model of specific interest to the present paper is the 
Hubbard model [55]. This lattice system is one of the 
most studied in research on strongly correlated systems, 
and is also the one which has received most attention 
in the context of static and time-dependent DFT stud- 
ies. For the non-equilibrium case, it is described by the 
Hamiltonian 



H= -t 



(RR')(j 



R 



Ra 



(8) 

Eq. ([8| is a direct generalization of the usual Hubbard 
Hamiltonian to the inhomogeneous, time-dependent case, 
and in presence of spin-dependent potentials |56| . Such a 
Hamiltonian offers one of the simplest descriptions of the 
competing behavior between the itinerant and localized 
behavior of electrons in the presence of interactions. In 
Eq. M, the term W{t) = Yl,Ra'^R<'('^)^R<y describes a 
local (in space and time), spin- and/or time-dependent 
potential (r denotes the time variable). For convenience, 
in W{t) we separate the static and time-dependent parts: 
WRa = e_R(T + Vf{cr{T)- lu the static case, all vr„(t) = 0; 
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furthermore, in spin-independent formulations, en^ and 
fijo-(T) are independent of u. Finally, the standard single- 
band homogeneous Hubbard Hamiltonian |55| is recov- 
ered when Ur = U and wjia ~ 0. 



A. General aspects of ground state lattice DFT 

The use of the Hubbard model in the context of static 
DFT was introduced in a comparative study of many- 
body and DFT Fermi surfaces (57j. However, a DFT for- 
mulation based on the local lattice occupation numbers 
nn had already been introduced by the same authors in 
earlier work |58l EHI , to study the XC discontinuity in a 
semiconductor model. Later on, a more general formula- 
tion of static lattice DFT was presented ISUj , and an anal- 
ysis of the local density approximation (LDA) was per- 
formed for the ID Hubbard model (where an exact solu- 
tion based on Bethe-Ansatz is possible [filj). At the same 
time, other formulations were proposed, which use the 
lattice one-particle density matrix 7^^/ as the basic vari- 
able [621 ES] ■ Further significant progress within static 
lattice DFT was made in Ref. [51]. In this work, a LDA 
based on the Bethe-Ansatz solution (henceforth denoted 
BALDA) for the v^c was proposed, suitable for a DFT 
of the inhomogeneous ID Hubbard model. An analytical 
paramatrization of the XC energy and potential was also 
provided. In a series of works [MUSS] , the BALDA for Vxc 
was tested and benchmarked against exact diagonaliza- 
tion, density-matrix renormalization group (DMRG) and 
quantum Monte Carlo calculations and shown to attain 
an accuracy of the order of a few percent for energies, 
particle densities and entropies. 

In this paper, we do not consider magnetic effects, and 
limit ourselves to a discussion of a spin-independent DFT 
and TDDFT for the Hubbard model. In static DFT, and 
in standard notation, we can write for the ground-state 
total energy [571 110] : 



Ey[n] = To[n\ + S//[n] + E^dn] + ^ext(»)"i, 

i 



(9) 



where Ugi:* is the static external field (in the notation of 
Eq. (|8|, Vextii) = In Eq. Q, Ui = J^a'^i'^^ '^tiile 
To[n] and Eh — | X^i ^i^l ^''^s, respectively, the non- 
interacting kinetic energy and the Hartree energy. To 
perform a local density approximation, E^c is obtained 
from a homogeneous Hubbard model, chosen as a suitable 
reference system: 



E 



H- 



(10) 



To obtain the exchange-correlation potential v^ci one per- 
forms the derivative of the XC energy/site Cxc = Exc/L 
with respect to the density (in the general case, a func- 
tional derivative should be considered): 



dexcjn, U) 
dn 



Due to electron-hole symmetry, 

exc{n,U)^exc{2~n,U), (12) 

in the entire density range [0, 2]; thus 

Vxc{n) = -Wxc(2 - n). (13) 

Finally, a local density approximation is performed: 

Vxcii) = Vxcirii). (14) 

To ease the numerics when Vxc is discontinuous (see be- 
low), one can use a slightly smoothened version of Vxc 
near n = 1. The v^c thus obtained can be used in ground- 
state DFT-LDA calculations, which amounts to solving 
self-consistently the Kohn-Sham (KS) equations 



(15) 



where t denotes the matrix for the single-particle hop- 
pings among nearest-neighbor sites, and ip^ is the K-th 
single-particle KS orbital, with m = J2Keocc\fi^i^)\'^ ■ 
The effective potential matrix is diagonal: (yKs)ii = 
VKsii) = vnii) +Vxcii) + Vextii), with vnii) = ^UiU^ 
being the Hartree potential. 



B. DFT for the Hubbard model and dimensionality 

ID Hubbard model.- At half filling, i.e. when n = 1, 
for the infinite, homogeneous Hubbard model, the exact 
ground-state energy (obtained by the Bethe-Ansatz) is 

m- 

rn 2/3 . TT -4Jo(a:)Ji(x) 
^^^'^^ = 'V""(^)-X' ^[l-fe.M^./2)] - 

The interpolation formula proposed for the XC en- 
ergy/site at general densities is [M] 

, 2/?(C/) . , TTTi , 4 . .Tin. Un^ 

exc{n,U) = sm — — ) -I- - sm — ) — , 

TT P[U) n 2 4 



(17) 

where /? depends on U but not on the density, and is 
determined by Eq. (16 1. It is easily seen that f3{U = 
0) = 2, and /3{U — ?> 00] = 1. This gives for Vxc- 



,,n<l 



(n)--2cos( — )-f2cos( — )-— (18) 



and v"^^{n) = —Vxc{'^ — n). With the parametrization 
of Eq. ( 18 ) for Vxc, the discontinuity of Vxc at half-filling 
is Avxc = 4cos(7r//3(J7)) -I- U f65|. This expression re- 
covers correctly the exact asymptotic limits A'^^*''^ = 
{8/n)VUe-^''^^ and A^-*°° ^ U - 4 + 81n(2)/C/. For 
U > 2 the parametrization of Eq.(18) agrees within few 



(11) 



percent with the full numerical Bethe-Ansatz solution. 
For U < 2, Eq. (18) is not equally accurate. We fi- 



nally mention that, very recently, a parametrization of 
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the XC potential has become available also for the spin- 
dependent case |69l. 

2D Hubbard model.- In 2D, the Hubbard model has been 
investigated via DFT on the graphene lattice. In this 
case, the expression for the XC energy/site is [7D] 

e,e(n, U) = 0(6-"^' - l)e-(^l"-i|-'^)' . (19) 

The parameters a, b, c, d were determined by a fitting 
procedure to the ground-state energy calculations, 
performed with the exact diagonalization method. 

3D Hubbard model.- To date, the only case considered 
in the literature is that of the simple cubic lattice [7T] . 
where the ground-state energy of the uniform system was 
computed within dynamical mean field theory (DMFT) 
[721 [73] . The short summary below follows rather closely 
the one originally given in 7X • 

Initially devised for lattice models with infinite connec- 
tivity (where the self-energy S is local in space), DMFT 
deals non-perturbatively with correlations; in 3D, it can 
be seen as an approximation scheme with a local self- 
energy. DMFT maps a Hubbard model on a simple cu- 
bic lattice onto a local problem representing one of the 
lattice sites (site 0) hybridized with a bath (the rest of 
the lattice) [73] . A Hamiltonian description of the local 
problem can be recovered, by identifying the site with 
the impurity site of an Anderson impurity model (AIM), 
and introducing auxiliary degrees of freedom {ajj^i- 

y-AiM = y-imp + ^ [efeoL "'k:<^ + ^ofc (oLcoo- + li-c.) , 

(20) 

where T-Limp = UnQ^riQi — /i(not + '^o^)- The value of 
the parameters {efc,Vofe} in the auxiliary Hamiltonian 
T~Laim are determined self-consistently, i.e. when the 
impurity single-particle Green's function G(zw„) in Mat- 
subara space [74] becomes identical to the local lattice 
Green's function with identical self-energy E(ia;„): 

G(iw„) = j deD{e)[ioJn -I- /i - e - S(icj„)]-\ (21) 

Here, S(iaj„) — Go^iii-^Jn) ~ G^^iii^n) is the local Dyson 

equation, (iujn) = jw„ + ^^ - Y.k iuj„-ek ' ^^"^ ^('^) 
is the non-interacting lattice density of states. Once at 
self-consistency, the relevant quantities {e.g., the double 
occupancy (no^no^), the total energy, etc.) are finally 
computed. It is perhaps worth mentioning that schemes 
other than DMFT could be used to determine the Cxc of 
the homogenous system, e.g. the Quantum Monte Carlo 
method or the Gutzwiller approximation (such work is 
currently under way). 

In Fig. [l]we present results for e^^^^'^ and v^^'^^^ 
for several U values. On increasing U , the curvature 
of Cxc changes, and for n « 1 a cusp develops above a 
critical value C/^^°'*. As a consequence, a discontinuity 
shows up in Vxc- This latter feature is the manifestation 



U=8 U=12 U=24 U=36 




FIG. 1: (Color online) Exchange-correlation energies e^ic (top) 
and potentials v^c (bottom) for the homogeneous 3D Hubbard 
model, as a function of the density, for several values of the 
interaction U . 



of the Mott- Hubbard metal-insulator transition within a 
DFT description. This behavior is rather different from 
what is observed in the ID Hubbard model, where v^c 
is discontinuous for any C/ > [M]. The XC potentials 
shown in Fig. [l]will be utilized in Sect. [V] 



C. Lattice TDDFT 

We have seen that DFT can be a viable route to the de- 
scription of the ground-state properties of strongly corre- 
lated systems such as the Hubbard model. A subsequent 
step would be to look at the dynamical properties of the 
Hubbard model via TDDFT. Similarly to the case of real- 
istic materials, this could be done at two different levels: 

1) One can invoke a linear response treatment to exter- 
nal fields, using the apparatus of many-body perturba- 
tion theory, starting from the Kohn-Sham Hamiltonian; 

2) In presence of strong time-dependent fields, one could 
instead resort to a full dynamical description based on 
TDDFT, by propagating in time the KS equations. 

For ID Hubbard-type Hamiltonians, initial work in the 
first direction was made in Refs. [75^77] . In these pa- 
pers, the main emphasis was on the study of the Hubbard 
model as a simplified system to gain insight into the gen- 
eral aspects of TDDFT, specifically the kernel /^c, and 
no real-time dynamics was performed. TDDFT in the 
linear response regime was also applied to a ID spinless 
fermion model with nearest-neighbour interactions [78], 
where comparisons between exact results and an LDA 
based on the Bethe-Ansatz were reported. Less is known 
for D > 1, and we will present below new results for 
the linear response of the 3D Hubbard model within the 
TDDFT-ALDA. 

A TDDFT approach to the real-time dynamics of the 
Hubbard model out of equilibrium was initially consid- 
ered in |39| . In this work, the exact many-body non- 
equilibrium dynamics of small ID Hubbard chains of dif- 
ferent length and particle density was performed, and 
from it, via a reverse-engineering procedure, the exact, 
time-dependent XC potential was obtained, by propa- 
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gating the time-dependent Kohn-Sham equations: 

{i+VKs{T))ip^{T)^idr(pn{T). (22) 



Eq. (22) is the generahzation of Eq. (15) to the time- 



dependent case. In [55], exact resuhs for the density and 
the XC potential were compared to those from an approx- 
imate Vxc obtained within the ALDA, and an analysis of 
non-local adiabatic effects in v^c was carried out. For 
the ALDA, in analogy to the static case, the ID Hub- 
bard model was chosen as the reference systems, so that 



.ALDA 



RALDA 



(23) 



where wf/^^'^(n) is given by Eq. (18). In Ref. 



the treatment was limited to spin-compensated systems, 
and the spin-dependent case was presented in |79| , where 
TDDFT results were compared to those obtained with 
the time-dependent DMRG method. 

An early rigorous discussion of the formal aspects of 
TDDFT can be found in [35], where an example of par- 
ticle density which is non wo-representable on the lat- 
tice was provided. In [381 it was shown that a tight- 
binding Hamiltonian may in general be unable to re- 
produce a specific given density. The case considered 
was that of a dimer with one orbital/site. The issue of 
wo-representability was further analyzed in [39l |80] . As 
shown in [80], for a ID lattice with L sites, the con- 
dition \Sk\ < 2y/nknk+i, with Sk = Z]i=i '^i (this ex- 
presses current conservation), is necessary and sufficient 
for WQ-representability for = 1 particle. In a lattice 
Hamiltonian, the spread of the eigenmodes scales with 
the hopping term t. According to the time-energy un- 
certainty relation, the range of possible response times 
scales also with \t\. Thus, the larger the hopping, the 
faster the system can respond in time, and the broader 
the range of WQ-representability in time |80j . The above 
condition for Sk is only a necessary one for A^ > 1 [39] . 
also when Hubbard interactions are taken into account. 
For the special case of a dimer, and using the Cauchy- 
Schwartz inequality, one can show that the exact density 
of a Hubbard dimer is WQ-representable |39j . 



D. Developments in lattice DFT and TDDFT, and 
applications. 

• Electric polarizahility.- Very recently, lattice DFT has 
been used to determine the polarizahility a of the ID 
Hubbard model [ST. Results for a from lattice DFT are 
in good agreement with DMRG ones. In particular, the 
response of the XC potential is in the same direction of 
the perturbing potential. Also, the possibility of dealing 
in lattice DFT with large samples made possible to ex- 
amine the scaling properties of a |81j . 

• Magnetic properties. - lattice DFT has also been used to 
study the Hubbard model on the graphene lattice, where 
it correctly describes the ground-state spin configuration 
of large graphene clusters (flakes) |70j . 



• Quantum Information and Entanglement entropy.- A 
key quantity to quantum information phenomena is the 
entanglement, or quantum correlations, of a system. En- 
tanglement has also been studied in the context of quan- 
tum phase transitions. A lattice DFT approach to the en- 
tanglement entropy of the Hubbard model was presented 
in [68 and, very recently, an explicit density functional 
for the entanglement was also provided [HI] . An TDDFT- 
ALDA approach to the time evolution of entanglement 
entropy was used in [83] , to discuss the expansion of ul- 
tracold fermion clouds in optical lattices. 

• Quantum transport.- A combination of DMRG and lat- 
tice DFT was used to gain insight into the exact ground- 
state XC functionals for a correlated-electron model sys- 
tem coupled to external reservoirs [Ml. The specific 
model considered was the so-called interacting resonant 
level model, for which DMRG and DFT conductances 
were found in good agreement. Very recently, a compar- 
ative assessment of lattice DFT and DMRG in transport 
has been provided in [85] . while an application of the 
spin-dependent BALDA |69j to quantum transport can 
be found in [5S] . 

The role of a discontinuity in a time-dependent descrip- 
tion of quantum transport has been recently examined 
within lattice TDDFT- ALDA [87]. Following the time- 
evolution of a single Anderson impurity attached to two 
biased leads, a dynamical notion of the Coulomb blockade 
was provided. Accordingly, Coulomb blockade manifests 
itself as a periodic sequence of charging and discharging 
of the nanostructure. This emerges also from a descrip- 
tion based on TDCDFT [T7] . In passing, we mention that 
ground-state current DFT has been also applied to lat- 
tice models [88l|89]. For ID Hubbard rings threaded by 
a magnetic flux, the effects of lattice impurities on the 
persistent currents and on the Drude weights has also 
been examined [89]. 

• Connection to many-body approximations, and criti- 
cal analysis of the local density approximations.- The so- 
called variational approach to TDDFT gS] HH] has the 
advantage of a systematic inclusion of many-body con- 
tributions in the XC potential. In this way non-locality 
in space and memory effects can be properly included, 
once Vxc is retrieved from the many-body self-energy via 
the time-dependent Sham-Schliiter equation [90]. This 
well established relation between many-body perturba- 
tion theory and TDDFT on the Keldysh contour was 
numerically examined in [113| . by looking at exchange- 
correlation potentials obtained via time-dependent re- 
verse engineering using the time-dependent densities 
from the KBE. The performance of two self-interaction 
correction schemes for the ID Hubbard model has been 
scrutinized in [91], while a shortcoming of the BALDA 
was pointed out in [92j . specifically its capability to re- 
produce correctly the Friedel oscillations in the inhomo- 
geneous ID Hubbard model. Finally, the lattice DFT for 
the ID Hubbard model has been used to illustrate the 
Hohenberg-Kohn mapping between densities and ground- 
state wave functions in terms of the metric-space aspects 
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of the Hilbert space |93) . 

• Cold atoms and dynamical quenches.- In optical lat- 
tices [M], it is possible to study fermionic and bosonic 
atoms with repulsive and/or attractive interactions and 
(because an accurate tunability of the lattice parame- 
ters is possible) often more directly and easily than in 
solid-state experiments. Trapped ultracold atoms on an 
optical lattice permit to study different ground-state sce- 
narios for the Hubbard model. Lattice DFT and lattice 
TDDFT have been used to investigate these systems. For 
example, an extensive study of the ground-state proper- 
ties of trapped repulsive fermions in a ID Hubbard lat- 
tice was provided in |95j. Here, a comparison of lattice 
DFT and quantum Monte Carlo density profiles was per- 
formed, and an detailed microscopic picture of the conse- 
quences of the interplay between particle-particle inter- 
actions and confinement was given. 

A similar analysis within lattice DFT has been car- 
ried out in [S^ for the case of a ID Hubbard model with 
attractive interactions. For lattice TDDFT, a first ap- 
plication to fermions on optical lattices was to compare 
(spin-dependent) ALDA and DMRG dynamics after the 
quench of a localized external perturbation |79j . A more 
recent study of the dynamics after a local quench can 
be found in ^7\. The effect of a global quench in ID, 
i.e. the removal of a parabolic trapping potential, was 
studied in |83j . where quenches of different speed (from 
adiabatic to sudden) were applied, showing a dynamical 
self-stabilization of the Mott insulator phase in the den- 
sity profile. In this work, the role of entanglement as a 
pointer of dynamical phase changes was also analyzed, 
and the system's thermalization was discussed from a 
TDDFT perspective. Finally, the beating and damping 
regimes of the Bloch oscillations in a 3D Hubbard optical 
lattice has also been investigated with TDDFT \7T\ . 



E. TD[C]DFT on the lattice: Runge-Gross 
theorem and connections to lattice TDDFT 

There is at present no formulation of the Runge-Gross 
(RG) theorem of TDDFT for the lattice. The reason for 
this has been explicitly discussed in [T7] and, in essence, 
it has to do with the fact that a reductio ad ahsurdum 
strategy, along the lines of the RG theorem, is not pos- 
sible. This is because, in the proof for the lattice, terms 
related to the one-particle density matrix appear, which 
do not carry a definite sign. This observation corrob- 
orates the fact that, for lattice Hamiltonians, one can 
actually find examples which contradict the uniqueness 
of the correspondence between densities and potentials. 

If one adopts the bond current as the basic variable, 
a one-to-one mapping can be established [15 betweeen 
current densities and Peierls phases [551 US] of the bond- 
hopping terms. In this way a rigorous formulation of 
TDCDFT on the lattice becomes possible pMTT] . either 
in analogy with the original RG theorem |15j . or via a 
reformulation based on non-linear differential equations 



methods [T^. 

In addition to establish a rigorous connection also for 
the densities, a lattice TDCDFT formulation has the fur- 
ther advantage of being valid in the presence of magnetic 
fields of a general kind, which enter via the Peierls phases 
(more precisely, the Peierls phases describe the effect of 
the vector potential A). This is not the case of TDDFT; 
it can be easily shown that, by a local gauge transfor- 
mation, the scalar onsite potentials of TDDFT can be 
removed from the Hamiltonian while, at the same time, 
introducing Peierls phases ipij — ai — aj in the hopping 
term between sites i and j. The sum of phases of such 
specific form over a lattice plaquette adds up to zero, 
and corresponds to a null line integral of A over a closed 
loop, i.e. to a zero magnetic flux across the plaquette. 
For these cases, i.e. in the absence of magnetic fields, 
TDDFT could be used; otherwise, TDCDFT should be 
considered. 

In conclusion, care must certainly be exerted in using 
lattice TDDFT when dealing with issues for which the 
uniqueness of the density-potential mapping is relevant. 
It should, however, still be possible to use lattice TDDFT 
on heuristic grounds, even if, at present, the practical 
implications of the lack of a rigorous formulation largely 
remain to be seen. 



IV. KADANOFF-BAYM DYNAMICS 

The main aim of the non-equilibrium Green's function 
(NEC) technique gQl EH [Ml [M] is to obtain expecta- 
tion values of single-particle operators and the total en- 
ergy of an interacting system subject to an external field. 
The NEC technique is used in a wide variety of fields, 
treating both real and model systems |102Hlllj . The 
method has the advantage of having a direct connection 
to MBPT, in which conserving [ID] approximations of 
increasing complexity can be constructed systematically 
and where memory effects |100| llOlj are automatically 
built in. However, when the NEC technique is used in 
conjunction with MBPT, the method has also a few lim- 
itations, some of which are practical in character. Before 
discussing further this point, we proceed now to a brief 
summary of the approach. 

A. Formalism 

The central object in the NEG echnique is the path- 
ordered, one-particle Green's function, 

Gil,2)^-^{T^[i,il)^\2)]), (24) 

where, in more detail, 1 denote single-particle space/spin 
and time labels, riaiti. Here, orders the times ^1,^2 
on the Keldysh contour [41] 7 {ti , t2 can be either real or 
imaginary) , and the field operators are in the Heisenberg 
picture; the brackets () denote averaging over the initial 
state (ground state or thermal equilibrium). 
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The G is determined by the Kadanoff-Baym equations, 
its equations of motion. Specializing to time ti, we have 

m 

iidt,-h{ti))Gih,t2) = Sil2)+ I ^{h,t)G{tM)dt. 



(25) 

Here h is the single-particle Hamiltonian and S, the ker- 
nel of the integral equation, is the self-energy, treated 
within a given many-body approximation (MBA). Eq. 
( 25 1 and the one corresponding to t2 are solved numeri- 



cally, and for non-isolated systems, E contains an addi- 
tional contribution, the embedding self-energy Eemb, to 
treat the coupling system-environment |103| 11091 1113j . 

The MBA:s we consider in this paper are the second 
Born, the GW and the T-matrix approximations (BA, 
GWA and TMA respectively). The self-energy in the 
BA includes all terms up to second order. The GWA 
[115| amounts to add up all the bubble diagrams which 
give rise to the screened interaction, W = U + UPW, 
where P(12) = G'(12)G(21). In this case the self-energy 
is I](12) = G'(12)W^(12). In a spin-dependent treat- 
ment of the TMA [TT^ ITTB] one constructs the T by 
adding up all the ladder diagrams, T = $ — $J7T, where 
$(12) = G(12)G(12). The expression of the self-energy 
then becomes E(12) = / J7(13)G(43)r(34)C/(42)d34 (for 
further details see e.g. |113| ). 

All these approximations are conserving [?D|, i.e. the 
macroscopic conservation laws (for the energy, the num- 
ber of particles, etc.) hold, which is a key requirement 
for time dynamics. The fact that the approximations are 
conserving does, however, not guarantee the soundness of 
other important features, as shown in the next Section. 



B. Some drawbacks of KBE+MBPT 

On the practical side, a first limitation is that the 
Green's function is a two-point function and thus scales 
quadratically with the propagation time. A second one 
is that the KBE must be solved self-consistently with a 
non-local kernel at every point in time. Currently, this 
is viable for very small systems, whilst larger ones can 
only be treated in the steady state regime (where only 
one time variable is required). 

On the more fundamental side, the NEG approach 
within MBPT manifests other undesirable traits, for ex- 
ample perturbation schemes may break down for strong 
interactions. However, even for converging schemes, un- 
physical features can emerge, which can sometimes play 
an important role. For example, in MBPT, the ground- 
state spectral function of a finite system contains in- 
finitely many poles, whereas in the exact solution the 
number of poles in finite. Moreover, in the KBE-based 
time dynamics, the MBA:s introduce an artificial corre- 
lation induced damping which, in some cases, can com- 
pletely dominate the evolution. For finite systems, such 
damped dynamics results in a steady (i.e. never decay- 



ing) state. In extended systems, a similar behavior is 
observed, but with a difference, namely arbitrarily long- 
lived states occur which, however, eventually decay into 
a unique steady state. This behavior is certainly artificial 
in finite systems, while its physical soundness in extended 
systems is much harder to establish. 

Another drawback concerns KBE and correlation func- 
tions: It is well known that single-particle Green's func- 
tions give access to certain two-particle properties such 
as, e.g., the total energy. In particular, they can be 
used to determine the density-density correlation func- 
tion at equal times, e.g. the so-called double occupancy, 
(fifi^nfi^) |1171I118] . However, within conserving MBA:s, 
correlation functions may violate important properties 
such as positiveness |117j . 

In spite of these possible drawbacks, the KBE, together 
with a MBPT approach to the self-energy, offer the great 
advantage of incorporating non-locality and memory ef- 
fects on equal footing. For this reason they will be used 
in the next Section to benchmark our TDDFT results. 



V. SCRUTINIZING THE LINEAR RESPONSE 
FORMALISM AND THE ALDA IN LATTICE 
MODELS 

This section presents original results for three different 
topics. The first one is an illustration, by means of a 
simple model system, of the problems one may encounter 
when using the LSS within the linear response formalism 
of TDDFT. 

We then present two applications of lattice TDDFT, 
which offer a clear illustration of the limitations (but also 
of the positive aspects) of the ALDA. In the first case, 
we study the density response function for the 3D Hub- 
bard model, and we use an f^J^^^ obtained from DMFT 
(Sect 



IIIBl. The functional form of the kernel /, 



ALDA 
xc 



is also discussed. In the second example, we study the 
non-linear, TDDFT- ALDA response of a finite cluster to 
strong, time-dependent perturbations. The Vxc used in 
the ALDA is also in this case obtained from DMFT. The 
TDDFT results will be compared to exact ones, and to 
those from the Kadanoff-Baym dynamics introduced in 
Sec. im 



A. Inversion of the non-interacting response 

It was long hoped that the exact-exchange approxi- 
mation within TDDFT would provide a reasonable XC 
kernel which is both non-local and energy dependent and 
thus constituting a significant step beyond the usual adi- 
abatic approximations. As mentioned in Sect. |Tlj these 
hopes have recently evaporated, the main reason be- 
ing zeros in the non-interacting response function. This 
problem has actually been pointed some time ago [54j . 
In the spirit of the model investigations at the heart of 
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the present work we will here demonstrate the problem 
in some model cases. 

The non-interacting response function xo can always 
be written as 



with a vanishing eigenvalue we find that 



Xo(i',r';w) 



2c.J,(r)/*(r') 
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(26) 



where the index g is a double index q — {k,k') and 
— Ck' — e/c is an excitation energy given by the en- 
ergy difference between an unoccupied (cfe/) and an oc- 
cupied (e/c) eigenvalue of the basic one-electron Hamilto- 
nian. The functions fq are excitation functions consisting 
of the corresponding products between an unoccupied, 
(pl,{r), and occupied, ipk{^), orbital of the same Hamil- 
tonian. Due to the orthogonality of the orbitals {(fik}, 
the functions fq always integrate to zero, resulting in the 
well known property of xo of a finite system 



Xo(r,r')dV = 0. 



(27) 



This means that the response function xo always has a 
zero eigenvalue at all frequencies, which is a reflection of 
the physical fact that there can be no charge response to 
a constant potential. Unfortunately, the response xo also 
has other zero eigenvalues at specific frequencies. If we 
first consider a two-electron problem, the lowest orbital 
ipo is doubly occupied and all other orbitals are unoc- 
cupied. The complete set of orbitals ipk are, of course, 
linearly independent, a fact which is not altered by re- 
moving the lowest member ipo from the set. Neither is 
this fact altered by multiplying all the remaining orbitals 
in the set by the same function tpQ. Consequently, in this 
case, all the excitation functions are linearly independent. 
Assuming for a moment that xo tias an eigenvector u(r) 
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FIG. 2: (Color online) Inverse of the eigenvalues of xo, for 
a six-site, spin-compensated chain with 2 and 4 electrons, as 
function of the frequency. The longer vertical ticks on the 
frequency axis denote the excitation energies (see main text). 
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(28) 



where the coefficients Uq are the projections of the eigen- 
function u onto the excitation functions fq. But, since 
the excitation functions are linearly independent, all coef- 
ficients in this linear combination of the fq:s must vanish 
and, since the frequency dependent factors are non-zero, 
we are led to the conclusion that all the coefficients Uq 
must vanish. But, again, the only function u which is 
orthogonal to all the excitation functions fq is a constant 
and we have recovered the already known case. Thus, 
in the two-electron case, there can be no additional zero 
eigenvalues of xo ■ 

Proceeding now to more than two electrons, we have 
several sets of excitation functions, one for each occupied 
orbital. Although all functions within a particular set are 
linearly independent among themselves, we find it highly 
unlikely that the conjunction of all sets are linearly in- 
dependent. Indeed, every set is almost a complete set by 
itself. One would therefore expect to be able to expand 
a member of one set in the combined functions of two 
other sets. Consequently, by varying the frequency and 
the coefficients Uq in the vanishing sum above, Eq.(28l, 
one would expect to be able to arrive at a particular 
vanishing linear combination of the excitation functions 
fq. In the spirit of the present paper, we will illustrate 
these points on a simple non-interacting linear chain of 
atoms with one orbital per site, no one-site energies, and 
only nearest neighbour hopping- the same for each spin 
channel. 

Choosing four atoms in the chain and two electrons 
with opposite spin, the relevant Hamiltonian is 4x4 and 
we have only three excitation functions with correspond- 
ing excitation energies - from the ground state to each of 
the three unoccupied states. As it turns out and is easily 
verified on the back of an envelope, the three excitation 
functions in a four-dimensional space are definitely lin- 
early independent and there are thus no additional zero 
eigenvalues of xo ■ 

Turning then to four electrons, two for each spin chan- 
nel, we obtain four excitation functions from each of the 
two occupied states to each of the remaining two unoc- 
cupied states. As easily verified, two excitation functions 
are identical and so are their corresponding excitation 
energies. In the sum above, Eq.(28 these two functions 



will contribute equally and both are included by multi- 
plying one of the terms by two. Any linear chain has, 
of course, mirror symmetry around the mid point, a fact 
which substantially reduces the necessary algebra - espe- 
cially for the few-site cases. The response function xo 
splits into a sum of an even and an odd contribution 
and so do all the excitation functions. In the present 
case, the two even excitation functions are identical and, 
consequently, the even part of x consits of only one term 
which cannot be orthogonal to anything but the constant 
vector. The two odd excitation functions are easily seen 
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local: 



FIG. 3: (Color online) The ALDA exchange-correlation kernel 
fx<t'^^ for [/ = 8 and U = 2A obtained by differentiating v^c 
with respect to the density. 



to be linearly independent by inspection. Thus, also in 
the four-electron case, there are no extra zero eigenvalues 
of xo as long as there are four sites. 

Including six sites in the model, we immediately ap- 
proach the above discussed general behaviour of xo- 
Looking first at the two-electron case, we obtain five ex- 
citation functions and energies. The calculations now 
become a little too complicated to carry out by hand but 
are still rather trivial. The six-by-six determinant of the 
five excitation functions plus the constant vector turns 
out to be non-singular, thus demonstrating all excita- 
tion functions to be linearly independent. Hence, there 
are no additional vanishing eigenvalues of xo, in keeping 
with the general discussion above. 

Having instead four electrons on the six sites gives us 
eight excitation functions of which two turn out to be 
identical. But there is no way the remaining seven can 
be linearly independent in a six-dimensional space. As 
discussed above, we would then expect that a particular 
choice of frequency could result in coefficients giving rise 
to a vanishing linear combination of excitation functions 
and thus a zero eigenvalue of xo- In Fig- |2]we have plot- 
ted the inverse of the five non-trivial frequency dependent 
eigenvalues of xo for the even and the odd channels - for 
both the two-electron case and the four-electron case. 
We see that these inverse eigenvalues are all nice smooth 
functions of the frequency in the former case. We also see 
that two eigenvalues in the four-electron case pass zero at 
particular frequencies away from the excitation energies 
- one in the odd channel and one in the even channel. 
These results corroborate the correctness of the general 
discussion above. 



B. Linear response in the ALDA for the 3D 
Hubbard model. 

In analogy to the continuum case, the general expres- 
sion for the XC kernel fxc on the lattice is 



(29) 



In the ALDA, an ordinary derivative in the density is 
performed, and the space and time dependence is entirely 



/,1,^^^(R,R';i,t') 



_ dvxc[n-R,{t)) 

-<,(nR(0)(SRR'<5(t-t'), (30) 



and fxc = fxc^°^{q,uj), i.e. fxc only depends of the 
(uniform) ground-state density of the system. In the 
Hubbard model of Eq. (|8|, the bare interaction U is 
also local in time and space. Accordingly, for the 3D ho- 
mogeneous Hubbard model, and in the (q, w) space, the 
ALDA density-density response function x{l,^) is 



xo(q,w) 



1 - (t^ + /a;c)xo(q,w) 



(31) 



Here, x(q)W) is the Fourier transform of the retarded 
response function 



x(R,i) = -i^(i)(["R(i),"o(0)]) 



gs, 



(32) 



and h^it) = nu{t) - {nn.it)), with hn{t) = nRt(t) + 
ftRiit)- Furthermore, xo is the response function of the 
Kohn-Sham system, 



xo(q,w) 



(27r)3 



^3j„ ^f{£u) - »F(ek+q) 
Ck - Sk+q + 0J + ir]' 



(33) 



where ek = —2t{coskx + cos ky + cos kz) is the single 
particle energy dispersion in the 3D simple cubic lattice, 
np{e) is the Fermi function (we work at zero tempera- 
ture) and the integral in k is performed over the first 
Brillouin zone. 

Before presenting the random-phase-approximation 
(RPA) and TDDFT-ALDA results for x, it is useful to 
discuss briefly the features of fxc in the ALDA. In Fig. 
§ we show /^/^^ as a function of the density n. Since 
fxc^^i'Ti) = fxc^^^C^ - n)i we can consider results for 
n < \. Furthermore, the results fxc for n < 0.2 and 
n > 0.95 exhibit significant noise, and thus are not dis- 
played. In Fig. |3j we see clearly one of the interesting 
features of fxc^ , namely the XC kernel can be either 
negative or positive [119] . This is especially evident for 
U = &, and can also immediately be gathered from the 
results for Vxc in Fig. Jl] Thus we see that, depending on 
the band filling, f^c^ can reduce or reinforce the effect 
of the bare interaction U in Eq. (31 1 . This is at variance 



with the usual case, where the XC kernel tend always 
to induce an effective interaction which is reduced with 
respect to the bare one. The other important aspect in 
the XC kernel is its behavior at n = 1. Looking again 
at Fig. n] we note that the discontinuity in v^^^^'^ for 
U > (in our case U = 24) will introduce a (Dirac's 

delta- like) spike in the XC kernel at n = 1. This singu- 
lar behavior would be absent for U < J7*^°" in 3D, but 
always present in ID where Vxc is discontinuous for any 
value of U. 

Results for x(q, i^) for three values of q are shown in 
the three panels of Fig. [4] The systems density that we 
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FIG. 4: (Color online) The imaginary part of — xCq, ^) for 
U = 8, as a function of a;, at three high symmetry points in 
the Brillouin zone. The chosen q values are X = (7r,0,0), 
M = (tt, tt, 0) and R = (vr, n, n) (note the different vertical 
scale for the three panels). The average density is n — 0.5, 
i.e. quarter filling. In each panel, the solid curves refer to 
the Kohn-Sham, RPA and ALDA cases (red, blue and green, 
respectively). The effect of f^^c for this filling is to give a 
reduced effective interaction, with the ALDA peak shifted to 
the left of the RPA one. For the point X we also show results 
for another filling (n = 0.85, dashed curves), where fxc is 
positive. In this case, the ALDA peak (green curve) is at the 
highest energy. 



consider is quarter filling, i.e. n = 0.5, and U = 8. In 
each panel, there are three solid curves which represent 
the imaginary part of the Kohn-Sham (red curve) , RPA 
(blue curve) and TDDFT-ALDA versions of -x(q,w) 
(the dashed curves in the top panel will be discussed 
momentarily). For the RPA and ALDA curves, we see 
sharp structures, outside the xo continuum region, cor- 
responding to the zeros of 1 — ([/ + fxc)xo{c[,uj)- Such 
peaks describe the plasmonic features of our system. We 
mention in passing that we have performed several test 
of our numerics, including a (successful) verification the 
of the f-sum rule, 

^^(q) = -;^y Im(xo(q,w))'^rft^ = 

(34) 

It interesting to note that, for the solid curves, the RPA 
peaks are always at higher energy than the ALDA ones. 
This is easily explained considering that at quarter filling, 
and for U = 8, f^c^^^ is negative (Fig. [s]), thus screening 
the bare U. However, it also clear that for densities close 
to half-filling, say n = 0.85, the XC kernel at U — 8 
is positive, giving ALDA peaks at energies higher than 
the RPA ones. This situation is shown in the top panel 
of Fig. |4] for one value of q, and it corresponds to the 
dashed curves in that panel. 



A similar behavior can be observed in real space. In 
Fig. [5j which refers to the case of quarter filling and U = 
8, the four panels show the dependence of Imx(R, cj) on 
the distance |R| along the x-axis (as one goes away from 
R — 0, the response function diminishes quite rapidly). 
The principal effect of the correlations is to move the 
structures in Imx to higher energies with respect of the 
Kohn-Sham results. Furthermore, we can notice how the 
ALDA profiles are always shifted at lower energies with 
respect to those of the RPA, due to the negative value of 
fxc^^^ when n = 0.5 and f7 = 8. 
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FIG. 5: (Color online) The linear response xi^!^) as a func- 
tion of iu for (7 = 8, for different values of R. 

In general, the pair correlation function (and thus the 
density response function) is a careful indicator of the 
accuracy of a many-body approximation. Thus, look- 
ing at properties which involve x can provide a severe 
test to assess the performance of the ALDA. For exam- 
ple, the density-density response function can be used 
to compute the ground-state energy of the system. More 
specifically, via the connection to the dynamical structure 
factor S{q,uj), one can determine the local double occu- 
pancy — {fi-B.-]-n-R.i) ■ We have used x in the RPA and 
the ALDA to compute ^r, and the results are reported 
m Table in Two general features can be noted. The first 
is that in some cases the RPA and ALDA results for 
d-R, are negative. Since dji is manifestly positive, this is 
an adverse outcome of the ALDA, and is rather general 
in character. Indeed, it is a long established fact that, 
in the RPA |120| . pair correlations for the electron gas 
at metallic densities are strongly negative at short dis- 
tances |121) : more recently, the same problem has been 
also experienced within the self-consistent GWA and BA 

The second general trait in Table|l]is related to the sign 
of fxc'- at quarter filling, where the fxc is negative, we get 
an effective lower interaction, which tends to increase the 
ALDA double occupancy with respect to the one from the 
RPA. On the other hand, at n ^ 0.85, the fxc for [/ = 8 
is positive, and thus instead increases the interaction, 
which results in a lower double occupancy. 
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Parameters 


doMFT 


dKS 


d-RPA 


dALDA 


[/ = 8,n = 0.5 


0.036 


0.062 


-0.010 


-0.007 


[7 = 24,n = 0.5 


0.016 


0.062 


-0.047 


-0.035 


U = 8,11^ 0.85 


0.114 


0.178 


+0.072 


+0.059 



TABLE I; Double occupancies d — {nfn_i_) obtained from 
DMFT, compared against results from linear response. Re- 
sults from xo, Xrpa and from xalda are shown. 

As a conclusive remark, while our results for the double 
occupancy point out a shortcoming of the ALDA, they 
also indicate that the problem with the sign of d-R could 
be less acute than in the continuum case. In fact, at 
least for the cases we considered, double occupancies in 
the ALDA are only slightly negative and, most likely, this 
is a consequence of the natural cut-off introduced by the 
lattice at short distances. 
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FIG. 7: (Color online) Time-dependent density, no(r), of the 
central (interacting) site in a 5x5x5 simple cubic cluster sub- 

ject to a Gaussian potential Wg(r) = — 5e 2 (thin black 
dashed). The curves correspond to exact (thick solid orange), 
ALDA (thin solid black), BA (red dashed), GWA (blue dot- 
dashed) and TMA (solid green). The number of particles per 
spin of the equivalent cluster is 2. In a) f/ = 8 and in b) 
U = 24. 



C. Exact, ALDA and Kadanoff-Baym real-time 
dynamics in a small cluster. 

The system we choose to compare TDDFT-ALDA, ex- 
act, and Kadanoff-Baym results is a simple cubic clus- 
ter with 5"^ sites and open boundary conditions (Fig. 
|6] left). We choose the cluster to be highly inhomo- 
geneous, by having a single interacting impurity in the 
center R = 0, e.g., Ur = USuo in Eq. We also 

set w_R,(r) = wo{T)Siio. Due to cubic symmetry, only 
Nsym = 10 out of 125 one-particle eigenstates, those with 
non-zero amplitude at i? = 0, determine the static and 
time-dependent density at i? = 0, making the size of 
the exact configuration space manageable |124j . in the 
guise of an effective 10-site cluster (Fig. [6] right). Fur- 
thermore, if Efl^^o = 0, further use of symmetry gives 
Nsym = 7. In this case, the shape of the cluster corre- 
sponds to the R = impurity directly connected with 
to all the other 6 sites. However the R ^ Q states are 
not connected directly to each other, but only to the 
impurity. In the ground state, the way Ng electrons 
in the cluster distribute between the Nsym active and 
125 — Nsym spectator states corresponds to the exact 
many-body eigenstate with lowest energy. As a overall 




5x5x5 site cluster effective cluster 



FIG. 6: (Color online). Original 125-site cluster (left) and 
its effective image (right). In the original cluster, all sites are 
related by symmetry to the 10 sites explicitly shown. Accord- 
ingly, the site is the same in the two clusters. 



remark, we note that the simplification due to symme- 
try is very peculiar to our choice of a local interaction 
at a single site, and permits us to obtain exact dynami- 
cal results in 3D for a quite large cluster. Furthermore, 
such simplification remains in the presence of a time- 
dependent perturbation localized at the impurity site. 

In Fig. [7|we study the time-dependent density of the 
interacting impurity when subject to a slowly varying 
pulse for different values of the interaction strength. This 
is a situation in which the ALDA is expected to be ap- 
propriate. In the ground state, before the perturbation 
has been introduced, the systems are in the low density 
regime. The reason of this choice is that of the three 
MBA:s we consider, one of them, the TMA, performs 
rather well at these densities |111| . while the BA and 
GWA are not good in any of the regimes we considered. 

In Fig. [7j where the interaction strength is weak 
and the external field is slowly varying, we see that the 
ALDA, as well as the KBE+MBPT, give an overaff good 
description, both in the time window where the pertur- 
bation is actually present, and afterwards, when the sys- 
tem is performing free oscillations. When the interac- 
tion strength is geared up in Fig. [Tj^ we clearly see that 
the ALDA is superior to all the KBE results from the 
MBA:s, while the TMA performs better than the BA and 
the GWA. It is fair to say that overall ALDA provides a 
good description for these slow external fields. 

In Fig. |8j we show the exact, ALDA, and TMA results 
for case of a moderate interaction strength and a step per- 
turbation (BA and GWA, not presented, perform quite 
poorly). In this case, the shortcomings of the ALDA are 
evident. For example, the ALDA does not reproduce the 
long-time oscillations properly: in Fig. |8^, the ALDA 
performs well for r < 4, but after that, ALDA suffers an 
increasing dephasing with respect to the exact and TMA 
solutions. On the other hand, the TMA and the exact 
curve are in very good mutual agreement for this situa- 
tion, (see e.g. the region r w 14), which corresponds to 
a highly non-adiabatic, but weak, perturbation. 
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FIG. 8: (Color online) Time-dependent density, no(r), of the 
central (interacting) site in a 5x5x5 simple cubic cluster sub- 
ject to a step potential, Wg{T) = WoO{t). The curves corre- 
spond to exact (thick orange), ALDA (thin black) and TMA 
(green). The number of particles per spin of the equivalent 
cluster is 1. In a) Wo — —0.2 and in b) Wo — —2. 



For a stronger step perturbation (Fig. 8p), the 
KBE-^MBFT approach, represented here by the TMA, 
is able to describe the history dependence quite well, and 
significantly better than the ALDA, which becomes un- 
reliable. However, when both the interaction is strong 
and the external field is highly non-adiabatic, both the 
ALDA and the KBE-^MBPT fail to describe the tran- 
sient dynamics (not shown). 

To summarize, the ALDA, in suitable conditions, i.e. 
for slow perturbations, can perform well. However, we 
wish to remark that the ALDA we have considered here 
accounts for the effect of correlations in a non pertur- 
bative way, since it is obtained from DMFT reference 
calculations. Otherwise, for systems with much stronger 
interactions, a failure of ALDA should be expected also 
for external fields which vary slowly in time. 



VI. CONCLUSIONS 

The application of time-dependent density-functional 
theories to the non-equilibrium dynamics of lattice mod- 
els has recently received some attention in the scientific 
community. Conceptually firm ground has now been es- 
tablished for approaches which use the current density 
(or, more properly, the lattice bond current) as basic 
variable. By contrast, the early TDDFT-like approaches 
for lattice systems made use of the density as the basic 
variable. On the lattice, a density-based approach lacks 
a completely rigorous formulation, but it appears possi- 
ble to consider it on heuristic grounds, in the sense that 
the missing rigor should not be of significant practical 
consequence. 

It is in this spirit that in this work we have adopted 
lattice TDDFT as an investigative tool to discuss two 
general aspects of TDDFT, which are present also in the 
continuum formulation. These are the linear response 
formalism in TDDFT, and the adiabatic local density 
approximations (ALDA) . 

Instrumental to this strategy, the first part of our pa- 
per has been devoted to short reviews of linear response 



theory, lattice TDDFT, and also of the non-equilibrium 
Green's function method: the latter method was used to 
benchmark the results from lattice TDDFT. 

In the second part of this contribution, where our new 
results are presented, we have studied the linear and non- 
linear response of ID and 3D lattice Hubbard-type mod- 
els. We find that the limitations encountered in the con- 
tinuum case by the TDDFT linear response persist in 
simple ID lattice model systems. Furthermore, an ALDA 
treatment for the linear response of the 3D homogenous 
Hubbard model exhibits an important pitfall which is 
also present in the continuum case, namely short-range 
correlations are incorrectly accounted for by the ALDA, 
albeit at a lesser extent than for the electron gas. 

We have also shown that the sign of the ALDA XC 
kernel can change depending on the values of U and the 
density, and that for large enough values of U, fxc can 
exhibit a singular behavior at half-filling (due to the dis- 
continuity in Vxc^ which in turn is indicative of the Mott- 
Hubbard metal-insulator transition). Such characteris- 
tics of fxc are clearly visible in the energy dependence 
of the density response function. To our knowledge, the 
qualitative features of f^,^^^ just mentioned have no 
correspondence in the available XC kernels for contin- 
uum case. On general grounds, we expect they should be 
robust against the inclusion of non-local, non-adiabatic 
effects. 

There have been attempts to go beyond the ALDA 
and simultaneously include non-locality and frequency 
dependence in the kernel f^c- Several of these attempts 
have originated in the LSS or, equivalently, in the varia- 
tional approach starting from the Klein functional, and 
have, therefore, led to severe difficulties as discussed in 
the introduction. The problem is associated with use of 
the LSS or, within the variational appoach, with starting 
from the Klein functional. The "perpetrator" is the dou- 
ble inversion of the xo within these formulations. We can 
already anticipate that this particular problem will not 
be present if one instead choses to start from the more 
sophisticated Luttinger and Ward (LW) functional |125j . 
Whether or not the LW functional will give rise to other 
problems remains to be seen. Work along these lines is 
underway. 

In the non-linear case, our results can be summarized 
as follows. An ALDA which accounts for the effect of 
correlations in a non-perturbative way, can provide a 
satisfactory description of the non-equilibrium dynam- 
ics induced by slow perturbations. When non-local and 
non-adiabatic effects play an important role, i.e. for fast 
perturbations, the ALDA is clearly insufficient. This is 
especially evident if one compares the exact, TDDFT- 
ALDA and Kadanoff-Baym dynamics. In the latter, ef- 
fects beyond the ALDA are present, and are seen to be 
necessary to get a good agreement with the exact data. 

In conclusion, both in the linear and non linear cases, 
our results for lattice models confirm that going beyond 
the ALDA is certainly necessary in many instances, and 
that lattice models can be a useful way to scrutinize fun- 
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damental open issues in TDDFT and to explore possible 
avenues for improvement. 
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